/**
 * @file : RHS.h
 * @author: fenghuwuyu
 * @version 1.0
 * @mail: 3148107006@qq.com
 * @date: 2021年04月20日 星期二 21时15分42秒
 * @brief:
 */

#ifndef RHS_H
#define RHS_H
#include <iostream>
#include <vector>
#include <cmath>

class ProductRHS
{
	public:
		ProductRHS(){};
		~ProductRHS(){};
		virtual std::vector<double> RHS(std::vector<double> &_U) = 0;
};


class ThreeBody_RHS:public ProductRHS
{
	public:
		ThreeBody_RHS(){};
		~ThreeBody_RHS(){};
		std::vector<double> RHS(std::vector<double> &_U);
	private:
		const double _mu = 1 / 81.45;
};

std::vector<double> ThreeBody_RHS::RHS(std::vector<double> &_U)
{
    std::vector<double> du(6);
    du[0] = _U[3];
    du[1] = _U[4];
    du[2] = _U[5];
    double t2 = _U[0] + _mu - 1.0;
    double t3 = _U[1] * _U[1] + _U[2] * _U[2] + t2 * t2;
    double t4 = sqrt(t3 * t3 * t3);
    double t5 = _U[1] * _U[1] + _U[2] * _U[2] + (t2 + 1.0) * (t2 + 1.0);
    double t6 = sqrt(t5 * t5 * t5);
    du[3] = 2.0 * _U[4] + _U[0] - _mu * t2 / t4 - (1.0 - _mu) * (t2 + 1.0) / t6;
    du[4] = -2.0 * _U[3] + _U[1] - _mu * _U[1] / t4 - (1.0 - _mu) * _U[1] / t6;
    du[5] = - _mu * _U[2] / t4 - (1.0 - _mu) * _U[2] / t6;
    return du;
}

#else
// DO NOTHING
#endif

